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Most of the liquid-state theories, including glass-transition theories, are constructed on the basis 
of two-body density correlations. However, we have recently shown that many-body correlations, in 
particular bond orientational correlations, play a key role in both the glass transition and the crys- 
tallization transition. Here we show, with numerical simulations of supercooled polydisperse hard 
spheres systems, that the lengthscale associated with any two-point spatial correlation function 
does not increase toward the glass transition. A growing lengthscale is instead revealed by consid- 
ering many-body correlation functions, such as correlators of orientational order, which follows the 
lengthscale of the dynamic heterogeneities. Despite the growing of crystal-like bond orientational 
order, we reveal that the stability against crystallization with increasing polydispersity is due to an 
increasing population of icosahedral arrangements of particles. Our results suggest that, for this 
type of systems, many-body correlations are a manifestation of the link between the vitrification 
and the crystallization phenomena. Whether a system is vitrified or crystallized can be controlled 
by the degree of frustration against crystallization, polydispersity in this case. 



I. INTRODUCTION 

Amorphous materials have been of prime importance 
in our technology for millennia, from antique glass works 
to fashionable phones made of metallic glass. One of 
the new frontiers of the amorphous technology is in the 
design of amorphous drugs [1, 2], better absorbed by our 
metabolism with less side effects, that would be stable 
at room temperature. The main obstacle is a lack of our 
basic understanding of the physics of the glass transition, 
without any operative consensus despite half a century 
of intensive research [3, 4]. 

When cooled below its freezing temperature while 
avoiding crystallization, a liquid becomes supercooled. 
Upon further cooling, the dynamics slows down by many 
orders of magnitude leading to a material that is mechan- 
ically a solid without long range positional order, thus 
called amorphous. It is now known that the dynamics in 
a supercooled liquid is heterogeneous, with a length scale 
that grows when approaching the glass transition [5, 6] 
(see [4] for a review). The lengthscale defined by the 
dynamical heterogeneity is not static (one-time spatial 
correlation) but dynamic (two-time spatial correlation). 

The existence of a static (structural) length that would 
grow and accompany the dynamic heterogeneities is still 
not clear in the general case. However, in a class of sys- 
tem that includes polydisperse hard spheres, we have 
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shown [7] that some medium range bond orientational 
order reminiscent of the crystal exists in the supercooled 
liquid and grows toward the glass transition in the same 
way as the dynamical heterogeneity. The presence in 
glassy materials of structures locally reminiscent of crys- 
tals has been confirmed recently in amorphous silicon [8] 
and in a metallic glass [9]. While bond orientational or- 
der is a member of the class of many-body correlations 
between neighboring particles, it is yet unclear if a sim- 
ilar lengthscale can be extracted from two-body correla- 
tion functions. This question is particularly important 
considering that the Mode-Coupling Theory (MCT) of 
the glass transition takes as input only two-body quanti- 
ties, and similarly modern spin-glass-type theories of the 
structural glass transition [10-12] are not taking explic- 
itly into account many-body correlations. 

At polydispersities over 6 ~ 7 %, the system needs to 
fractionate to crystallize [13]. What is the bond order of 
the reference crystal is then unclear and may challenge 
our scenario. This is a situation reminiscent of binary 
hard sphere systems of size ratio close to one [14, 15] 
where growing bond order have not been reported [16]. 
However, it is known that even in binary systems locally 
favored structures play a role in the slowing down of the 
dynamics in some cases [17, 18]. Interestingly, recent 
studies by Mosayebi et al. [19, 20] demonstrate that there 
are static growing lengths in binary Lennard- Jones and 
soft sphere mixtures. They estimated the static length 
by looking at the spatial correlation of the degree of non- 
affine deformation of inherent structures under shear de- 
formation. 

In the present study we will use polydisperse hard 
sphere systems, where we know how to extract mean- 
ingful many-body correlations, and look for a two-body 
quantity that would show the same behavior as the bond 
order. We will show that the two-body part of the free 
energy (which for hard potentials corresponds to the two- 
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body part of the structural entropy) is unable to capture 
medium range bond ordered regions or to yield corre- 
lation lengths meaningful from the point of view of the 
glass transition. We thus confirm the medium range crys- 
talline order scenario and test its robustness against in- 
creasing polydispersity. Since bond orientational order 
is directly linked to the underlying crystalline structures, 
we will then address the important question of what is 
the mechanism responsible for the avoidance of crystal- 
lization. We will show that in the metastable fluid phase 
crystalline packings are in competition with icosahedral 
packings, and that polydispersity acts in favor of the lat- 
ter ones. 

The paper is organized as follows. In Sec. II we present 
the details of the simulations and the order parameters 
considered in this work. In Sec. Ill the results are orga- 
nized into a study of the order parameters distribution 
(Sec. Ill A) and their static lengthscales (Sec. IIIB), and 
a method to determine the competition between crys- 
talline structures and icosahedral packings (Sec. IIIC). 
In Sec. IV we discuss our results. In Sec. V we summa- 
rize our work. 



II. METHODS 

A. Simulation method 

We run isothermal-isobaric NpT Monte Carlo simu- 
lations of N = 4000 polydisperse hard spheres. The 
diameters (a) follow a Gaussian distribution P(cr) = 
exp [— (a — <to) 2 /2 s 2 ] /V^rs, with polydispersity index 
A = s/cjq. In the following we fix the unit of length 
as do = 1 and the unit of energy so that the Boltzmann 
constant is unity, ks = 1. 



B. Estimation of two-point quantities: pair entropy 
and pair free energy 

Our aim is to compare the behavior of both two- 
point quantities and many-body quantities with increas- 
ing pressure. Due to the hard-sphere interaction, entropy 
is the only contribution to the system free energy. All 
two-pair correlation quantities are thus derived from the 
two-body excess entropy [21, 22], defined as 

s 2 = -^Jdr [g(r) log( 5 (r)) - g(r) + 1] . (1) 

In principle, S2 can be calculated separately for each 
particle i in the system. In practice, this requires time 
averages to compute the pair correlation function of each 
particle, #i(r), where the particle distribution around a 
particle i is averaged over short-time scales (/? processes). 
In Refs. [7, 23, 24] we averaged on times comparable or 
longer than the a relaxation, leading to a quantity that 
was trivially a reflection of the dynamical heterogeneity. 



Here we instead construct an approximate but instan- 
taneous S2(i) using the pair correlation function g(r) 

S2 (i) = ]T bfoi) lQ g(5(^)) " 9(nj) + 1] • (2) 
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This quantity is in very good agreement with the one 
obtained by calculating the radial distribution function 
for each particle, ^(r), averaged over times comparable 
to the f3 relaxation time. 

More rigorously, one can compute directly the free- 
energy of each configuration by measuring the free vol- 
ume of the particles, defined as the volume (v(i)) in which 
each sphere can freely move while holding all the other 
spheres fixed. It has been shown [25] that this free vol- 
ume is simply related to the pair free-energy (^2) by the 
following relation 

/2 = £/2(i) = -fc B T£log(t;(i)/A), (3) 

i i 

where A is the thermal de Broglie wavelength. To com- 
pute the free volume v(i) we follow previous studies [25]: 
first the Voronoi-diagram for each configuration is com- 
puted, and the polyhedron surrounding each particle is 
determined. To account for polydispersity we employ 
the radical Voronoi tessellation. The free volume of par- 
ticle i is then computed by shifting normally all the faces 
of the corresponding polyhedron by cr(i)/2 toward par- 
ticle z, and computing the new volume. In this way the 
volume v(i) represents the volume in which the excluded 
volume of particle i can move without leaving its Voronoi 
cell. This procedure is conducted independently for each 
particle and for each configuration. 

C. Estimation of many-point quantities: bond 
orientational order parameter analysis 

To study many-body correlations we use the lo- 
cal bond-order analysis introduced by Steinhardt [26], 
first applied to study crystal nucleation by Frenkel 
and co-workers [27]. The i-fold symmetry of a neigh- 
borhood around each particle i is characterized by a 
(2£ + 1) dimensional complex vector (qj) as q£ m (i) = 
NbfJ) ^2f=i^ ^m(ry), where i is a free integer parameter, 
and m is an integer that runs from m = —£ to m = £. The 
functions Y^ m are the spherical harmonics and f ij is the 
vector from particle i to particle j. The sum goes over all 
neighboring particles Nb(i) of particle i. Usually Nb(i) 
is defined by all particles within a cutoff distance, but 
in an inhomogeneous system the cutoff distance would 
have to change according to the local density. Instead 
we sort neighbors according to their distance from parti- 
cle z, and fix Nb(i) = 12 which is the number of nearest 
neighbors in icosahedra and close packed crystals (like 
HCP and FCC) which are known to be the only relevant 
crystalline structures for hard spheres. 
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In the analysis, one uses the rotational invariants de- 
fined as: 



\ 2Z + 1 ^ '^ m|2 ' 

\ m=-£ 



(4) 
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where the term in brackets in Eq. (5) is the Wigner 3-j 
symbol. In particular both crystalline and icosahedral 
neighborhood have high q§ (strong 6-fold symmetry), 
with the highest values for the latter. To detect specif- 
ically icosahedral order one prefers Wq, whose minimum 
value is obtained only by a perfect icosahedron. 

The identification of crystalline particles follows the 
usual procedure [27]. A particle is identified as crystal if 
its orientational order is coherent (in symmetry and in 
orientation) with that of its neighbors. The scalar prod- 
uct (q 6 (i)/|q 6 (i)|) • (qeCO/lqeCOl) quantifies this simi- 
larity. If it exceeds 0.7 between two neighbors, they are 
deemed connected. We then identify a particle as crys- 
talline if it is connected with at least 7 neighbors [27]. 
In a more continuous way, summing the contribution of 
all the bonds of a given particle, we define the "crys- 
tallinity" [28] 



N b (i) 



C(0 = £ 



qeW • qe(i) 
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(6) 



Alternatively, one can coarse-grain over the neigh- 
bors [29] 

1 / m \ 

Qim{i) = N ^ ~ x \qem(i) + ^m(i) I , (7) 

to suppress the signal from locally incoherent orientations 
(icosahedral order) [30] and use the resulting invariant Qq 
as an indication of crystallinity, more precisely, the de- 
gree of crystal-like rotational symmetry. Alternatively, 
we note that the shortcomings of non-coarse-grained or- 
der parameters in the identification of crystallinity can 
be addressed by Minkowski tensors [31]. 

Here we briefly consider the physical meaning of 
crystal-like bond orientational order parameters. The 
crystallization transition is characterized by the symme- 
try breaking of both orientational and translational or- 
der. We note that both C and Qq are good measures 
of bond orientational order, whereas the density or other 
two-body quantities are measures of translational order. 
It was shown [28] that in hard spheres crystallization is 
driven by fluctuations in bond orientational order and 
not by density fluctuations. Crystals continuously form, 
grow and melt in regions of high bond orientational order, 
which then effectively act as precursors for the crystal- 
lization transition. So C and Qg, while not being direct 
indicators for the presence of crystals, rather measure 



the tendency to promote crystallization. In Sec. IIIB we 
are going to show that indeed the lengthscale associated 
with bond orientational order fluctuations increases with 
supercooling. Then in Sec. IIIC we are going to study 
the mechanism by which the crystallization transition is 
avoided. 



D. Estimation of the correlation lengths of various 
quantities 

Finally we explain how to evaluate the correlation 
length of various order parameters. The calculation can 
be carried out both in real space and in Fourier space. 
While formally containing the same information, the 
Fourier analysis has some practical advantages over the 
real space analysis. In real space, the correlation function 
of any order parameter is a oscillating and rapidly de- 
caying function of r. The correlation length is obtained 
by fitting the envelope of the correlation function with 
a Ornstein-Zernike expression. This expresion is only 
asymptotic, so the two first peaks at short r should be 
omitted, and it is also rapidly decaying, so that the sta- 
tistical noise strongly affects the quality of the fit at long 
r. The problem is less severe for order parameters having 
a tensorial nature, and correlation lengths of crystal-like 
bond orientational order have been easily measured in 
real space in previous studies [7, 28, 30, 32, 33]. For 
example, the tensorial order parameter Q^m effectively 
correlates 7 scalar fields, allowing a sevenfold reduction 
of the noise. But the real space analysis requires much 
longer time averages for two-body correlation functions, 
both ji and s<i, which have a purely scalar nature. This 
problem can be overcome by calculating the correlation 
functions in Fourier space. These function are not oscil- 
latory at small q where we can expect a Ornstein-Zernike 
form, and thus much easier to fit unambiguously. So, to 
keep both two-body and many-body correlation functions 
consistent, we calculate all correlations in Fourier space. 
We explain in the following a straightforward procedure 
to extract correlation length from Fourier space analysis. 

At a given time step, for any scalar order parameter 
field x increasing with order, we compute a structure fac- 
tor keeping only the 10% most ordered particles (more 
on this choice below). This condition defines a thresh- 
old x*. The ensemble average of the thresholds (x*) 
are indicated on Fig. 2. Formally we define a function 
= &[x(i) — x*], where Q(x) is the Heavisides step 
function, and a four-point structure factor 

S x (q) = N- 1 «fi(q)fi(-q)> - | <tt(q)} 2 |), (8) 
where fi(q) is the Fourier transform of cu(i): 

^(q) = ^^(z)exp(-zq-r i ). (9) 



This structure factor is then ensemble averaged (still 
noted S x (q) for concision) over w 10 4 configurations. The 
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FIG. 1: Simulated state points. The blue filled circles (with 
blue curve) represent the equation of state for polydispersity 
A = 7%. The red open squares (with red line) instead are 
simulation points at the same pressure (f3pcr 3 — 23) but at 
different polydispersities A: from low to high volume fraction 
they correspond to A = 7%, 9%, 11%, 13%, and 15% respec- 
tively. Similarly the black open diamonds (with black line) 
are at f3po 3 = 17 and A = 0%, 4%, and 7%. 

case of order parameters decreasing toward ordering (i.e. 
52, wq) is trivially obtained by changing the sign. 



III. RESULTS 

Figure 1 shows the equation of state for the simu- 
lated state points. In particular we consider three data- 
sets. The first one (blue filled circles in the figure) cor- 
responds to simulations at a constant polydispersity of 
A = 7%. For each state point we run 8 independent 
simulation runs and extract configurations for the cal- 
culation of correlation lengths. The second data set 
(red open squares in the figure) are instead isobaric 
simulation (at Ppcr 3 — 23) with increasing polydisper- 
sity, A = 7%, 9%, 11%, 13%, 15%. The third data set 
(black open diamonds in the figure) are also isobaric 
simulations at fipcr 3 = 17 with increasing polydisper- 
sity, A = 0%,4%,7%. The two last data sets are used 
to study the mechanism by which crystallization is sup- 
pressed upon an increase of polydispersity, unveiling the 
role played by icosahedral arrangement of particles. 



A. Order parameter distribution and mobility 

We study systematically two-body (s2, /b) an d many- 
body (#6, W6> Q6: C) scalar order parameter fields for 
our highest pressure (ftpcr 3 = 25). For each of these pa- 
rameters one can define if a particle is "ordered" or not. 
Very negative values of S2 indicate low two-body configu- 
rational entropy and thus a priori some kind of ordering 
or stability. Similarly, high values of fa indicate high 
two-body free energy (low free volume). High values of 
Qq or C indicate locally crystalline environment (locally 
similar to HCP or FCC crystal), whereas very negative 
values of wq correspond to icosahedral packings. High 
values of q§ can indicate ambiguously either crystalline 
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FIG. 2: Probability distributions (red dashed line) and mobil- 
ity (blue continuous line) function of various order parameters 
at (3pa 3 = 25 and for a time difference corresponding to the 
a-relaxation. Top row: two-body excess entropy S2 (left) and 
pair free energy fa (right). Central row: local six- fold orien- 
tational orders q$ (left) and wq (right). Bottom row: coarse 
grained Q6 (left) and crystallinity C (right). Mobility is in 
unit of mean-square displacement. The shaded area shows 
the contribution from 10% of the particles having highest (for 
fa, <?6, Qq and C) or lowest (for S2, wq) value of the order 
parameter. 



or icosahedral environments. 

In Fig. 2 we show the "ordered" side of each parame- 
ter as a shaded area, and the probability distribution of 
this parameter as a red dashed line. Note that for any 
parameter, the maximum probability do not correspond 
to the ordered side. However the probability distribution 
decays more slowly on the ordered side, indicating the 
presence of rare but very ordered particles. 

By plotting the probability distribution function of the 
metastable fluid in the (Qg, fa) and (wq, fa) planes, Fig. 3 
shows the absence of linear correlations between fa and 
both Q 6 and wq. Since Qq identifies regions of high 
crystal-like bond orientational order and wq locates icosa- 
hedral arrangements of particles, it is clear that high fa 
regions are not associated with any of these structures. 
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FIG. 3: Correlation between two-body and many-body pa- 
rameters. Probability distribution functions in the /2-Q6 map 
(left) and in the f2-we map (right) for a metastable fluid with 
polydispersity A = 7% and pressure f3pa 3 = 23. There is no 
linear correlation between the parameters, meaning that high 
values of Qq (the crystalline regions) or low values of wq (the 
icosahedral regions) are not extremum values of fa. 



We checked in the same way that 52 is not correlated 
with many-body parameters. 

The absence of strong correlations between two-body 
quantities and both crystalline and icosahedral packings 
is also evident from the microscopic dynamics. To study 
the correlation between any scalar order parameter x and 
the displacement of the particles, we define the mobility 



Ar 2 (x = xo,t) 



•Z\\n{t)-ri(0)\\ 2 S{x{i)-x )< 
J26(x(i) - x ) 



(10) 

shown in Fig. 2 for a time difference corresponding to 
the a-relaxation time. Note that we use 5 functions of 
a constant finite width and thus the number of parti- 
cles involved in the average of Eq. (10) varies like the 
probability distribution, explaining the noise in the low 
probability regions. 

We found that for each parameter, its mobility de- 
creases with increasing order. The mobilities of bond- 
order quantities are flat in the disordered regions and de- 
crease when approaching the perfect structure (i.e. icosa- 
hedron for wq, crystal structures for Qq or C). By con- 
trast, the mobility of two-body order parameters tends 
to increase strongly in the disordered region and decrease 
less in the ordered region (it is almost fiat at high /2). We 
conclude that many-body quantities describe better the 
slowing down accompanying good local ordering, while 
two-body quantities are not clearly correlated to such 
slower structures. Note that in Fig. 2 both icosahedral 
packings (low wq particles) and bond-ordered crystalline 
regions (high Qq and C) are associated with slow dynam- 
ics. As was shown by some of us [30], the structures pri- 
marily responsible for the slowing down of the dynamics 
are the crystal-like particles, while icosahedral particles 
act to prevent the crystal nucleation process [28]. This 




FIG. 4: Structure factors collapse on the Ornstein-Zernike 
form for (top) fi and (bottom) Qq . Note that the fi structure 
factors are almost similar for all pressures considered. For 
Qq instead, the growth of the correlation length is evident 
already from the systematic change in the relative position of 
the nearest-neighbors peaks. 



will be shown in detail in Sec. IIIC. 



B. Length-scales 

As explained in Sec. II D, we estimate the correlation 
lengths of various quantities x in Fourier space not only 
for the two-body 82 and /2, but also for scalars derived 
from the multi-body bond orientational order, i.e. qs, 
wq, Qq and C, which allows us to have overall coherency 
of the length-scale analysis of all the parameters. 

Figure 4 shows the increase of S x (q) toward small 
wavenumbers (x = ji upper panel; x — Qq lower panel), 
which is well described by the asymptotic Ornstein- 
Zernike function in Fourier space: 



S x (q -»• 0) 



Xx 



1 + &Z 2 ' 



(11) 



where ^ x is the correlation length and Xx the susceptibil- 
ity of fluctuations of quantity x. In general, an indepen- 
dent determination of Xx is crucial for the fit [34]. How- 
ever here we deal with correlation lengths much smaller 
than the simulation box and both £ x and Xx can be re- 
liably estimated from a two parameter fit of S x (q). We 
note that the absence of the finite-size effects was con- 
firmed for this situation [7]. 

The dependence on the pressure of the resulting cor- 
relation lengths £ x is shown in Fig. 5. Most of the order 
parameters produce constant lengthscales, including not 
only the two-body quantities but also q$, Wq that are 
sensible to icosahedral order. The only growing lengths 
are extracted from measures of local crystal-like order, 
i.e. Qq and C. We confirm that the same results are ob- 
tained from real-space correlation functions (but in real 
space is possible to detect the growth also in q§, see [35]). 
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FIG. 5: Correlation length (£ x ) and susceptibilities (x x ) ex- 
tracted for two-body and many-body scalar order parameters, 
plotted as a function of the pressure. Only many-body corre- 
lation lengths are increasing (only slightly for £ WQ ) , while the 
lengthscale associated with the two-body quantities is almost 
constant. 



The absence of correlations for two-body quantities 
and the presence of a growing lengthscale for Qq and C 
are evident also by direct inspection of the particle con- 
figurations. Figure 6 plots the 10% most ordered par- 
ticles for the different order parameters at f3pa 3 = 23 
and A = 7%. The first row of Fig. 6 shows the absence 
of any appreciable correlation for both two-body quan- 
tities /2 and 82. The middle row shows that also the 
signal from icosahedral clusters (both q 6 and Wq have 
icosahedra as their extremum) display no appreciable 
correlation length, i.e. they form randomly and homo- 
geneously throughout the system. Only Qq and C (last 
row in Fig. 6) show clustering of the ordered particles on 
medium range. 

The length scales obtained from Qq or C in Fourier 
space and from the tensorial q6 or Q6 in real space (not 
shown) are similar and increase monotonically with pres- 
sure. Note that the scalar q§ (dominated by the icosa- 
hedral order) yields a very different correlation length. 
This is consistent with the spatial coherency (in orienta- 
tion) of crystal-like order that is missing in icosahedral 
order [7, 30]. Coherently, real-space correlation function 
(not shown) of fi (respectively s^) are perfectly identical 
at all pressures. 




Qe C 

FIG. 6: Visualization of the 10% most ordered particles de- 
fined by the various order parameters. All pictures correspond 
to a thin slice (5<j) of the same configuration at f3pcr 3 — 23 
and A = 7%. Only Qq and C show meaningful spatial fluctu- 
ations. We can also see anti-correlation between (Qq,C) and 
(q6,w 6 ). 



The choice of the threshold x* is a balance between 
taking in too many particles or too few. If too few (be- 
low 5%) S x is too noisy. If too many, the threshold does 
not discriminate between ordered and disordered parti- 
cles and S x tends to the trivial density S(q). We found 
that between 5% and 40% the absolute value of the length 
depends marginally on but its pressure dependence 
does not. We chose to use 10% across this paper because 
this value allows the easiest direct visualization on a sin- 
gle frame (Fig. 6). We checked the pressure independence 
of the two-body parameter's length with thresholds up to 
90%. 

The study of correlation lengths has shown that by 
increasing pressure, the range of crystal-like bond orien- 
tational order increases, driving the slowing down of the 
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FIG. 7: Effect of polydispersity on local structures at con- 
stant pressure (f3pcr 3 = 23): (top) detail of the distribution of 
wq (full distribution in inset) showing a small increase in the 
icosahedra population saturating around 10%; (bottom) dis- 
tribution of Qq showing a marked decrease in the crystallinity 
with polydispersity. 



system as shown in Refs. [7, 30]. Bond orientational or- 
der corresponds to orientationally ordered regions which 
spontaneously form in the metastable phase, fi is in- 
stead decoupled from the relevant structures involved in 
the transition, as was shown in Fig. 3. We now address 
the question of how the system avoids crystallization de- 
spite the growing lenghtscale of bond orientational order. 



C. Competition between icosahedral arrangements 
and crystal-like arrangements 

We will now focus on the state point at fipa 3 = 23 
at different polydispersities to study the mechanism by 
which polydispersity disfavors the crystallization transi- 
tion. 

In Fig. 7 we show the probability distributions for the 
order parameters Qq and wq at different polydispersities. 
It is immediately evident that, while bond orientational 
order is rapidly suppressed with increasing polydispersity 
(as shown in the suppressed signal at high Qg), particles 
in icosahedral environments are not disfavored by poly- 
dispersity. On the contrary the fraction of icosahedral 
particles increases with polydispersity, and saturates at 
around A = 10%. This is in agreement with recent evi- 
dence of increased icosahedral ordering with size dispar- 
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FIG. 8: Average crystallinity order parameter projected on 
the q4-q6 plane for the metastable fluid at f3pa 3 = 23 at A = 
7% (left) and A = 15% (right). Icosahedra appear in the top- 
left corner of each plot (in dark gray) and perfect FCC crystals 
would be in the top-right corner with C = 12 (in black). 
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FIG. 9: Polydispersity dependence of the correlation lengths 
at /3pa 3 = 23. The dominant crystalline length decreases, the 
icosahedral length increases, however they plateau well before 
crossing. The two-body length shows no indication of becom- 
ing dominant with increasing A, even decreasing slightly. 



ity in metallic glasses [36]. Intuitively, we can conclude 
that icosahedral order is more tolerant to size asymmetry 
(with the small particle usually sitting at the center of 
the icosahedral cage) than crystalline order is. 

Figure 8 shows that the metastable fluid distribution 
on the q^-qe plane, which is a convenient representation 
as perfect icosahedral packings sit on the top-left corner 
of the distribution, while perfect crystals on the top-right 
corner. The figure clearly shows that, by increasing poly- 
dispersity, the biggest change in the metastable fluid dis- 
tribution is the suppression of crystal-like regions (high 
values of q±, q§ and C, the black region), while icosahe- 
dral environments are slightly enhanced (high q$, low q^ 
and negative C, the dark gray region). 

The different effects of polydispersity on icosahedral 
and crystalline ordering are reflected in the different cor- 
relation lengths. Figure 9 shows the correlation length 
extracted from Q 6 , wq and fi as a function of poly- 
dispersity. While the correlation length for Qq (asso- 
ciated with crystal-like regions) decreases with increas- 
ing polydispersity, the one extracted from wq (associ- 
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FIG. 10: Average as a function of q§ for particles identified 
as "fluid" and "crystalline" according to the criteria outlined 
in Sec. IIC. The continuous red curve represents fluid particles 
in a system at f3pa 3 = 17 and A = 0%, while dashed blue 
curve represents fluid particles at the same pressure but at 
A = 4%. The thick black curve represents instead crystalline 
particles (this curve is less sensitive to polydispersity and it 
is reported once). Note that, while in the monodisperse case 
there is a interval in which the crystalline particles have 
higher orientational order, in the polydisperse case crystalline 
particles have lower orientational order at all <j). While the 
monodisperse case is easily crystallized in direct simulations, 
the polydisperse system always remains met ast able. 



ated with icosahedral regions) increases. However the 
two lengths are far from crossing and seem to saturate 
around A = 13 — 15%. The correlation length of ji is 
constant or even slightly decreasing with A, never tak- 
ing over the many-body lengths. It is thus clear that, in 
the present range of polydispersity, crystal-like bond or- 
der fluctuations is still the dominant contribution in the 
static (and dynamic [30]) properties of the system. 

In Ref. [28] we have shown that the competition be- 
tween crystalline packings and icosahedral packings can 
be studied via two-dimensional maps of translational vs 
orientational order. Orientational order is captured by 
#6? which is small for disordered arrangement of parti- 
cles and increases for both crystal-like and icosahedral 
particles. Translational order is instead measured with 
the local packing fraction, (/>, obtained by measuring the 
volume of the Voronoi diagram associated with each par- 
ticle. The calculation is straightforward: for each con- 
figuration, crystalline particles are identified with the 
method described in Sec. II C and the other particles are 
termed "fluid" . For each subset of particles, the average 
value of the volume fraction <p is calculated as a function 
of the order parameter q$, and plotted in Fig. 10. For 
each value of q& the map captures the average volume 
fraction </> of both crystalline and non-crystalline envi- 
ronments. 

In Fig. 10 we compare the curves at f3pa 3 = 17 and 
at different polydispersities (see also the black diamonds 
of Fig. 1): at 0% (monodisperse system, red continuous 
curve) and at 4% (blue dashed curve). The curve for the 
crystalline particles is similar at the two considered poly- 
dispersities and is reported once as the thick solid line. 
First we consider the monodisperse case. As shown in 
the figure, at low volume fraction, a particle in the fluid 



phase has higher orientational order than in the crys- 
talline phase. But at <p = 55.8% a crossover occurs and 
the crystal phase gains microscopic stability: a particle in 
a crystalline environment will have higher orientational 
order than a particle in the fluid phase at the same vol- 
ume fraction. This crossover marks the appearance in 
the fluid phase of the metastable crystals which continu- 
ously appear, grow and shrink, until eventually a crystal 
droplet reaches the critical size and starts the crystal- 
lization process. At a higher volume fraction (0 = 58%) 
another crossover occurs, with crystalline particles hav- 
ing less orientational order than particles in the "fluid" 
branch. These particles in the fluid phase at high g 6 
are easily identified as particles in icosahedral packings 
(they have a low value of wq). It is thus clear that icosa- 
hedral packings are competing with crystalline packings, 
eventually dominating at high 0. This scenario is con- 
firmed by looking at the polydisperse case (blue dashed 
line in Fig. 10). At polydispersity A = 4% the crys- 
talline branch always lays above the fluid one, meaning 
that crystalline environments, at any fixed 0, are not able 
to attain higher bond orientational order than other con- 
figurations. Microscopically the difference between the 
monodisperse and the polydisperse case is due to an in- 
creased population of icosahedral particles, which domi- 
nate the fluid branch for q& > 0.5. We have shown that, 
while for the monodisperse case there is a range of volume 
fractions where crystalline particles attain higher orienta- 
tional order than icosahedral arrangements, for the poly- 
disperse case icosahedral particles always attain higher 
orientational order. This is immediately reflected in di- 
rect simulations, as at this pressure it is not possible to 
crystallize simulations at A = 4%, while the monodis- 
perse simulations are easily crystallized [37, 38]. Since 
the diffusional dynamics of the system, which controls 
the kinetic factor of crystallization [39] , is approximately 
the same at different values of A [37], this difference in 
the crystallization behavior has to come from some struc- 
tural difference introduced by polydispersity, i.e. an in- 
creased population of icosahedral particles (see Figs. 7 
and 8). We also confirm that at polydispersity A = 7% 
icosahedral particles are favored over crystalline ones for 
all the pressures studied, in line with the observation of 
Fasolo and Sollich [40] that one-phase crystallization is 
suppressed at high polydispersity. 

We have thus provided direct evidence that icosahedral 
particles are responsible for the suppression of crystalliza- 
tion in polydisperse hard spheres. 



IV. DISCUSSION 

In the present study we have compared the evolution of 
static length scales in polydisperse hard spheres for both 
two-body and many-body correlation functions. While 
two-body correlation functions do not show any sign of an 
increasing length scale with pressure, we have confirmed 
that in the glass transition of polydisperse hard spheres 
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the relevant static length is the correlation length of the 
crystal-like bond order. We have also determined the rel- 
evant microscopic structures that are associated with the 
increasing lengthscale: they correspond to crystal-like en- 
vironments of particles and are characterized by slow dy- 
namics. We thus confirm that in the glass transition of 
polydisperse hard spheres the relevant static length is the 
correlation length of the crystal-like bond orientational 
order. 

The other relevant structure with slow dynamics are 
icosahedral packings of particles, but their lengthscale 
does not grow appreciably with increasing pressure. 
Icosahedral assemblies of particles are spatially uncorre- 
cted. While not having a direct role in the slowing down 
of the dynamics, we have shown that icosahedral packings 
are responsible for the avoidance of the crystallization 
transition with increasing polydispersity. The increase in 
polydispersity reduces the degree of crystal-like bond ori- 
entational order whereas enhances the icosahedral order 
(see Figs. 7 and 8). The former is crucial for triggering 
crystal nucleation [28], whereas the latter leads to the 
frustration against crystallization, a role that increases 
with polydispersity. None of these physical aspects of 
the system could be described by two-body quantities. 

We previously showed that spatial fluctuations of 
crystal-like bond orientational order are closely corre- 
lated with local dynamics: more ordered regions have 
slower dynamics [7, 30, 41-44]. Together with these 
results, we may say that it is many-body correlations, 
or crystal-like bond orientational ordering, that are the 
cause of slow dynamics and dynamic heterogeneity. This 
means that future theories of glass transition and crys- 
tallization should deal with many-body correlation ef- 
fects properly. The link between the symmetry of the 
relevant bond order parameter for describing structural 
ordering in a supercooled liquid and that for crystalliza- 
tion is another interesting point. This may be a direct 
consequence of the fact that glass transition is governed 
by the same free energy as that controlling crystalliza- 
tion [45-47]. This conjecture is further supported by the 
role of polydispersity in the glass-forming ability of hard 
spheres. 

The mechanism by which polydispersity increases the 
barrier for crystal nucleation may be two fold: (i) direct 
random disorder effect which destroys crystal-like bond 
orientational order in a supercooled liquid, which is the 
first step in crystal nucleation [28]; (ii) enhancement of 
icosahedral ordering with an increase in the polydisper- 
sity. It is known that size disparity between a parti- 
cle and its surrounding neighbors stabilizes icosahedral 
order [36]. Since the symmetry of icosahedral order is 
not consistent with that of the equilibrium crystal poly- 
morphs (fee and hep for this case), competing ordering 
toward these two mutually inconsistent symmetries leads 



to strong frustration effects on crystallization, as in the 
case of 2D spin liquids [43, 48]. The results shown in Fig. 
10 suggest that mechanism (ii) may be more relevant for 
the suppression of crystallization. 

Although we studied polydisperse hard spheres in this 
article, these frustration mechanisms should be relevant 
to many other glass-forming systems including metallic 
glass formers [9, 49]. 



V. SUMMARY 

In this article, we show firm evidence for the impor- 
tance of many-body correlations in glass transition phe- 
nomena for hard spheres liquids. This feature cannot be 
described by the standard liquid-state theories based on 
two-body correlation. This implies that, at high density, 
liquid state packing effects inevitably lead to many-body 
correlations, which play key roles in phenomena like the 
glass transition and crystallization. A physically natural 
order parameter to pick up these many-body correlations 
is the bond order parameter, whose importance has been 
well recognized and established for ordering transitions 
of hard disks in 2D [50]. We believe in the importance 
of incorporating many-body correlations into theories to 
describe both the glass transition and crystallization phe- 
nomena properly [45, 46]. 

Our study also indicates that there is an intrinsic 
link between crystallization and vitrification. Whether 
a polydisperse hard spheres system is crystallized or vit- 
rified can be controlled just by changing polydispersity, 
which affects extendable crystal-like bond orientational 
order and isolated icosahedral order in an opposite man- 
ner. We speculate that this direct link may exist for 
systems where crystallization does not involve phase sep- 
aration, in other words, as far as the two phenomena are 
described by the same free energy [45-47]. How univer- 
sal is this scenario needs to be checked carefully in the 
future. 
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